#######################################################
#####REPLICATION ARCHIVE FOR NUANCED ACCOUNTABILITY####
#########BRITISH JOURNAL OF POLITICAL SCIENCE##########
#######WARD DATA, SETUP & NATIONAL ELECTIONS ONLY######
###de Kadt & Lieberman, 2017

#######################################################
################ SETUP BEGINS HERE ####################
#######################################################
## USER: change working directory on line 56 ##########
###Packages
library(foreign)
library(stargazer)
library(Hmisc)
library(sandwich)
library(lmtest)

###Cluster Robust SEs###
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

###Correlation Table###
cor.prob <- function(X, dfr = nrow(X) - 2) {
  R <- cor(X)
  above <- row(R) < col(R)
  r2 <- R[above]^2
  Fstat <- r2 * dfr / (1 - r2)
  R[above] <- 1 - pf(Fstat, 1, dfr)
  
  cor.mat <- t(R)
  cor.mat[upper.tri(cor.mat)] <- NA
  cor.mat
}

###Working Directories###
## USER: Input your own working directory for files on the next line ##
dir = paste("####INPUT WORKING FILE DIRECTORY BETWEEN QUOTATION MARKS ####")

##

fig = paste(dir,"output/", sep="")
stgz =  paste(dir,"output/", sep="")

###Data Setup###
#Load dataset
dat = read.dta(paste(dir,"nuanced_accountability_ward.dta", sep=""))
#Set up covariates
cov = paste("+ unempfrac_ch10 + log_income_ch10 + log_pop_ch10 + femalefrac_ch10 + whitefrac_ch10 + employment_broad2001 + log_income2001 + log_pop2001 + sex2001  + white_frac2001 + share_area_all_tbvc_ta + opp_control")

#######################################################
################ OUTPUTS BEGIN HERE ###################
#######################################################

#####################
##### SUM STATS #####
#####################
dat.sub.sum = dat[which(dat$refuse_80th==1 & dat$water_80th ==1 & dat$flush_80th),
                  c("water_ch10", "flush_ch10", "refuse_ch10", "sd_means_ch10",
                    "lgeanc_vs_ch10", "npeanc_vs_ch10",  "log_income_ch10", "unempfrac_ch10", "log_pop_ch10", "femalefrac_ch10", "whitefrac_ch10", 
                    "log_income2001", "employment_broad2001","log_pop2001","sex2001","white_frac2001", 
                    "share_area_all_tbvc_ta")]
stargazer(dat.sub.sum, 
          title="Summary statistics ward-level aggregate data",out = paste(stgz,"sumstatsagg.tex",sep="\\"))

#####################
### CORRELATIONS ####
#####################
###Table 8 in the Appendix: Intercorrelations between service changes
dat.sub = dat[which(dat$refuse_80th==1 & dat$water_80th ==1 & dat$flush_80th),c("water_ch10", "flush_ch10", "refuse_ch10", "sd_means_ch10")]
dat.sub = na.omit(dat.sub)
cortab = cor(dat.sub)
rownames(cortab) = c("Water", "Toilets", "Refuse", "Mean SD")
colnames(cortab) = c("Water", "Toilets", "Refuse", "Mean SD")
filename = paste("intercorrelations.tex")
stargazer(cortab,  
          label = paste("tab:intercorrelations"),
          title = paste("Do changes in services intercorrelate?",sep=""),
          out = paste(stgz,filename,sep="\\"))

##############################
##### TEST FOR SELECTION #####
##############################
i = 1
out.mod = list() 
lm.out = list()
out.mod.comp = list() 
lm.out.comp = list()
for(iv in c("anc_vs1999", "anc_vs2000")){
  for(dv in c("water_ch10", "flush_ch10", "refuse_ch10", "sd_means_ch10")){
    dat.80th = dat[which(dat$water_80th==1 & dat$flush_80th==1 & dat$refuse_80th==1),]
    
      dat.80th.comp = dat.80th[which(dat$competitive==1),]

      dat.sub = dat.80th[,c(dv, iv, 
                            "employment_broad2001", "log_income2001", "log_pop2001", "sex2001", 
                            "white_frac2001", "share_area_all_tbvc_ta", "cat_b")]
      dat.sub = na.omit(dat.sub)
      
      dat.sub.comp = dat.80th.comp[,c(dv, iv, 
                                      "employment_broad2001", "log_income2001", "log_pop2001", "sex2001", 
                                      "white_frac2001", "share_area_all_tbvc_ta", "cat_b")]
      dat.sub.comp = na.omit(dat.sub.comp)
      
    
    lm.out[[i]] = lm(paste(dv, " ~ ", iv, " + employment_broad2001 + log_income2001 + log_pop2001 + sex2001 + white_frac2001 + share_area_all_tbvc_ta", sep=""), 
                     data=dat.sub)
    out.mod[[i]] = coeftest(lm.out[[i]],vcov = vcovCluster(lm.out[[i]], dat.sub$cat_b))
    lm.out.comp[[i]] = lm(paste(dv, " ~ ", iv, " + employment_broad2001 + log_income2001 + log_pop2001 + sex2001 + white_frac2001 + share_area_all_tbvc_ta", sep=""), 
                          data=dat.sub.comp)
    out.mod.comp[[i]] = coeftest(lm.out.comp[[i]],vcov = vcovCluster(lm.out.comp[[i]], dat.sub.comp$cat_b))
    i = i+1
  }
}
clse.mod = lapply(out.mod, `[`,,2)
clse.mod.comp = lapply(out.mod.comp, `[`,,2)

###Table 9 in the Appendix (combination of these two stargazer tables): Does baseline ANC vote share predict changes in SD?
filename = paste("endogenousselectionnoncomp.tex")
stargazer(lm.out[1],lm.out[5],lm.out[2],lm.out[6],lm.out[3],lm.out[7],lm.out[4],lm.out[8],
          se = c(clse.mod[1],clse.mod[5],clse.mod[2],clse.mod[6],clse.mod[3],clse.mod[7],clse.mod[4],clse.mod[8]),
          style = "ajps",
          column.separate = c(2,2,2,2),
          column.labels = c("Water", "Toilets", "Refuse", "Mean SD"),
          label = paste("tab:endogenous-selection"),
          #covariate.labels = covars,
          title = paste("Does baseline ANC vote share predict service change?",sep=""),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "table",
          out = paste(stgz,filename,sep="\\"))

filename = paste("endogenousselectioncompetitive.tex")
stargazer(lm.out.comp[1],lm.out.comp[5],lm.out.comp[2],lm.out.comp[6],lm.out.comp[3],lm.out.comp[7],lm.out.comp[4],lm.out.comp[8],
          se = c(clse.mod.comp[1],clse.mod.comp[5],clse.mod.comp[2],clse.mod.comp[6],clse.mod.comp[3],clse.mod.comp[7],clse.mod.comp[4],clse.mod.comp[8]),
          style = "ajps",
          column.separate = c(2,2,2,2),
          column.labels = c("Water", "Toilets", "Refuse", "Mean SD"),
          label = paste("tab:endogenous-selection"),
          #covariate.labels = covars,
          title = paste("Does baseline ANC vote share predict service change?",sep=""),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "table",
          out = paste(stgz,filename,sep="\\"))

#######################################################
################ MAIN RESULTS HERE ####################
#######################################################

###################################
###NATIONAL GOVERNMENT ELECTIONS###
###################################
i = 1
#General holders for results:
out.mod.full.nocov = list()
out.mod.full.cov = list()
out.mod.full.nocov.comp = list()
out.mod.full.cov.comp = list()
out.mod.full.interact.cov = list()
out.mod.full.interact.cov.comp = list()

out.mod.full.nocov.incumb = list()
out.mod.full.nocov.incumb.comp = list()
out.mod.full.cov.incumb = list()
out.mod.full.cov.incumb.comp = list()
out.mod.full.nocov.opp = list()
out.mod.full.nocov.opp.comp = list()
out.mod.full.cov.opp = list()
out.mod.full.cov.opp.comp = list()

out.mod.full.cov.nowc = list()
out.mod.full.cov.comp.nowc = list()
out.mod.full.cov.incumb.nowc = list()
out.mod.full.cov.incumb.comp.nowc = list()
out.mod.full.cov.opp.nowc = list()
out.mod.full.cov.opp.comp.nowc = list()

out.mod.full.cov.nokzn = list()
out.mod.full.cov.comp.nokzn = list()
out.mod.full.cov.incumb.nokzn = list()
out.mod.full.cov.incumb.comp.nokzn = list()
out.mod.full.cov.opp.nokzn = list()
out.mod.full.cov.opp.comp.nokzn = list()

out.mod.full.cov.noec = list()
out.mod.full.cov.comp.noec = list()
out.mod.full.cov.incumb.noec = list()
out.mod.full.cov.incumb.comp.noec = list()
out.mod.full.cov.opp.noec = list()
out.mod.full.cov.opp.comp.noec = list()

lm.out.full.nocov = list()
lm.out.full.cov = list()
lm.out.full.nocov.comp = list()
lm.out.full.cov.comp = list()
lm.out.full.interact.cov = list()
lm.out.full.interact.cov.comp = list()

lm.out.full.nocov.incumb = list()
lm.out.full.nocov.incumb.comp = list()
lm.out.full.cov.incumb = list()
lm.out.full.cov.incumb.comp = list()
lm.out.full.nocov.opp = list()
lm.out.full.nocov.opp.comp = list()
lm.out.full.cov.opp = list()
lm.out.full.cov.opp.comp = list()

lm.out.full.cov.nowc = list()
lm.out.full.cov.comp.nowc = list()
lm.out.full.cov.incumb.nowc = list()
lm.out.full.cov.incumb.comp.nowc = list()
lm.out.full.cov.opp.nowc = list()
lm.out.full.cov.opp.comp.nowc = list()

lm.out.full.cov.nokzn = list()
lm.out.full.cov.comp.nokzn = list()
lm.out.full.cov.incumb.nokzn = list()
lm.out.full.cov.incumb.comp.nokzn = list()
lm.out.full.cov.opp.nokzn = list()
lm.out.full.cov.opp.comp.nokzn = list()

lm.out.full.cov.noec = list()
lm.out.full.cov.comp.noec = list()
lm.out.full.cov.incumb.noec = list()
lm.out.full.cov.incumb.comp.noec = list()
lm.out.full.cov.opp.noec = list()
lm.out.full.cov.opp.comp.noec = list()

for(iv in c("water_ch10", "flush_ch10", "refuse_ch10", "sd_means_ch10")){ #Sets the four different Dependent Variables
  dat.80th = dat[which(dat$water_80th==1 & dat$flush_80th==1 & dat$refuse_80th==1),] #Subsets the data to exclude the top quintiles
  
  dat.sub = dat.80th[,c(iv, 
                        "npeanc_vs_ch10", "unempfrac_ch10", "log_income_ch10", "log_pop_ch10", "femalefrac_ch10", "blackfrac_ch10", 
                        "whitefrac_ch10", "indianfrac_ch10", "colouredfrac_ch10", "traddwell_ch10", "cat_b", "opp_control", "province2011",
                        "employment_broad2001","log_income2001","log_pop2001","sex2001","indian_frac2001","colored_frac2001","black_frac2001","white_frac2001",
                        "share_area_all_tbvc_ta","traddwell_proportion2001","refuse_perc2001","all_flush_perc2001","home_water_perc2001","sd_means2001",
                        "competitive")] #Subsets the data to only the relevant variables
  dat.sub = na.omit(dat.sub) #Applies na.omit
  
  dat.sub.comp = dat.sub[which(dat.sub$competitive==1),] #Creates subset that is only competitive
  
  dat.sub.comp = na.omit(dat.sub.comp)
  
  iv.int = paste(iv,"*opp_control", sep="") #Creates interaction term for analyses of heterogeneous effects
  
  if(iv=="water_ch10"){
    baseline = " + home_water_perc2001"
  } else if(iv=="flush_ch10"){
    baseline = " + all_flush_perc2001"
  } else if(iv=="refuse_ch10"){
    baseline = " + refuse_perc2001"
  } else {
    baseline = " + sd_means2001"
  }
  
  ##Main specifications
  lm.out.full.nocov[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, sep=""), data=dat.sub)
  lm.out.full.cov[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, cov, baseline, sep=""), data=dat.sub)
  
  lm.out.full.nocov.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, sep=""), data=dat.sub.comp)
  lm.out.full.cov.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, cov, baseline, sep=""), data=dat.sub.comp)
  
  lm.out.full.interact.cov[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub)
  lm.out.full.interact.cov.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp)
  
  ##By incumbency status
  lm.out.full.nocov.incumb[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, sep=""), data=dat.sub[which(dat.sub$opp_control==0),])
  lm.out.full.cov.incumb[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==0),])
  lm.out.full.nocov.incumb.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0),])
  lm.out.full.cov.incumb.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0),])
  
  lm.out.full.nocov.opp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, sep=""), data=dat.sub[which(dat.sub$opp_control==1),])
  lm.out.full.cov.opp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==1),])
  lm.out.full.nocov.opp.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1),])
  lm.out.full.cov.opp.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1),])
  
  ##Excluding particular provinces
  lm.out.full.cov.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$province2011!="Western Cape"),])
  lm.out.full.cov.comp.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$province2011!="Western Cape"),])
  lm.out.full.cov.incumb.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="Western Cape"),])
  lm.out.full.cov.incumb.comp.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="Western Cape"),])
  lm.out.full.cov.opp.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="Western Cape"),])
  lm.out.full.cov.opp.comp.nowc[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="Western Cape"),])
  
  lm.out.full.cov.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$province2011!="KwaZulu-Natal"),])
  lm.out.full.cov.comp.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$province2011!="KwaZulu-Natal"),])
  lm.out.full.cov.incumb.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="KwaZulu-Natal"),])
  lm.out.full.cov.incumb.comp.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="KwaZulu-Natal"),])
  lm.out.full.cov.opp.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="KwaZulu-Natal"),])
  lm.out.full.cov.opp.comp.nokzn[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="KwaZulu-Natal"),])
  
  lm.out.full.cov.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$province2011!="Eastern Cape"),])
  lm.out.full.cov.comp.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$province2011!="Eastern Cape"),])
  lm.out.full.cov.incumb.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="Eastern Cape"),])
  lm.out.full.cov.incumb.comp.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="Eastern Cape"),])
  lm.out.full.cov.opp.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="Eastern Cape"),])
  lm.out.full.cov.opp.comp.noec[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="Eastern Cape"),])
  
  #######CLUSTER STANDARD ERRORS:
  out.mod.full.nocov[[i]] = coeftest(lm.out.full.nocov[[i]],vcov = vcovCluster(lm.out.full.nocov[[i]], dat.sub$cat_b))
  out.mod.full.cov[[i]] = coeftest(lm.out.full.cov[[i]],vcov = vcovCluster(lm.out.full.cov[[i]], dat.sub$cat_b))
  
  out.mod.full.nocov.comp[[i]] = coeftest(lm.out.full.nocov.comp[[i]],vcov = vcovCluster(lm.out.full.nocov.comp[[i]], dat.sub.comp$cat_b))
  out.mod.full.cov.comp[[i]] = coeftest(lm.out.full.cov.comp[[i]],vcov = vcovCluster(lm.out.full.cov.comp[[i]], dat.sub.comp$cat_b))
  
  out.mod.full.interact.cov[[i]] = coeftest(lm.out.full.interact.cov[[i]],vcov = vcovCluster(lm.out.full.interact.cov[[i]], dat.sub$cat_b))
  out.mod.full.interact.cov.comp[[i]] = coeftest(lm.out.full.interact.cov.comp[[i]],vcov = vcovCluster(lm.out.full.interact.cov.comp[[i]], dat.sub.comp$cat_b))
  
  out.mod.full.nocov.incumb[[i]] = coeftest(lm.out.full.nocov.incumb[[i]],vcov = vcovCluster(lm.out.full.nocov.incumb[[i]], dat.sub[which(dat.sub$opp_control==0),]$cat_b))
  out.mod.full.cov.incumb[[i]] = coeftest(lm.out.full.cov.incumb[[i]],vcov = vcovCluster(lm.out.full.cov.incumb[[i]], dat.sub[which(dat.sub$opp_control==0),]$cat_b))
  out.mod.full.nocov.incumb.comp[[i]] = coeftest(lm.out.full.nocov.incumb.comp[[i]],vcov = vcovCluster(lm.out.full.nocov.incumb.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0),]$cat_b))
  out.mod.full.cov.incumb.comp[[i]] = coeftest(lm.out.full.cov.incumb.comp[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0),]$cat_b))
  
  out.mod.full.nocov.opp[[i]] = coeftest(lm.out.full.nocov.opp[[i]],vcov = vcovCluster(lm.out.full.nocov.opp[[i]], dat.sub[which(dat.sub$opp_control==1),]$cat_b))
  out.mod.full.cov.opp[[i]] = coeftest(lm.out.full.cov.opp[[i]],vcov = vcovCluster(lm.out.full.cov.opp[[i]], dat.sub[which(dat.sub$opp_control==1),]$cat_b))
  out.mod.full.nocov.opp.comp[[i]] = coeftest(lm.out.full.nocov.opp.comp[[i]],vcov = vcovCluster(lm.out.full.nocov.opp.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1),]$cat_b))
  out.mod.full.cov.opp.comp[[i]] = coeftest(lm.out.full.cov.opp.comp[[i]],vcov = vcovCluster(lm.out.full.cov.opp.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1),]$cat_b))
  
  out.mod.full.cov.nowc[[i]] = coeftest(lm.out.full.cov.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.nowc[[i]], dat.sub[which(dat.sub$province2011!="Western Cape"),]$cat_b))
  out.mod.full.cov.comp.nowc[[i]] = coeftest(lm.out.full.cov.comp.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.comp.nowc[[i]], dat.sub.comp[which(dat.sub.comp$province2011!="Western Cape"),]$cat_b))
  out.mod.full.cov.incumb.nowc[[i]] = coeftest(lm.out.full.cov.incumb.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.nowc[[i]], dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="Western Cape"),]$cat_b))
  out.mod.full.cov.incumb.comp.nowc[[i]] = coeftest(lm.out.full.cov.incumb.comp.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.comp.nowc[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="Western Cape"),]$cat_b))
  out.mod.full.cov.opp.nowc[[i]] = coeftest(lm.out.full.cov.opp.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.opp.nowc[[i]], dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="Western Cape"),]$cat_b))
  out.mod.full.cov.opp.comp.nowc[[i]] = coeftest(lm.out.full.cov.opp.comp.nowc[[i]],vcov = vcovCluster(lm.out.full.cov.opp.comp.nowc[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="Western Cape"),]$cat_b))
  
  out.mod.full.cov.nokzn[[i]] = coeftest(lm.out.full.cov.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.nokzn[[i]], dat.sub[which(dat.sub$province2011!="KwaZulu-Natal"),]$cat_b))
  out.mod.full.cov.comp.nokzn[[i]] = coeftest(lm.out.full.cov.comp.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.comp.nokzn[[i]], dat.sub.comp[which(dat.sub.comp$province2011!="KwaZulu-Natal"),]$cat_b))
  out.mod.full.cov.incumb.nokzn[[i]] = coeftest(lm.out.full.cov.incumb.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.nokzn[[i]], dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="KwaZulu-Natal"),]$cat_b))
  out.mod.full.cov.incumb.comp.nokzn[[i]] = coeftest(lm.out.full.cov.incumb.comp.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.comp.nokzn[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="KwaZulu-Natal"),]$cat_b))
  out.mod.full.cov.opp.nokzn[[i]] = coeftest(lm.out.full.cov.opp.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.opp.nokzn[[i]], dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="KwaZulu-Natal"),]$cat_b))
  out.mod.full.cov.opp.comp.nokzn[[i]] = coeftest(lm.out.full.cov.opp.comp.nokzn[[i]],vcov = vcovCluster(lm.out.full.cov.opp.comp.nokzn[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="KwaZulu-Natal"),]$cat_b))
  
  out.mod.full.cov.noec[[i]] = coeftest(lm.out.full.cov.noec[[i]],vcov = vcovCluster(lm.out.full.cov.noec[[i]], dat.sub[which(dat.sub$province2011!="Eastern Cape"),]$cat_b))
  out.mod.full.cov.comp.noec[[i]] = coeftest(lm.out.full.cov.comp.noec[[i]],vcov = vcovCluster(lm.out.full.cov.comp.noec[[i]], dat.sub.comp[which(dat.sub.comp$province2011!="Eastern Cape"),]$cat_b))
  out.mod.full.cov.incumb.noec[[i]] = coeftest(lm.out.full.cov.incumb.noec[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.noec[[i]], dat.sub[which(dat.sub$opp_control==0 & dat.sub$province2011!="Eastern Cape"),]$cat_b))
  out.mod.full.cov.incumb.comp.noec[[i]] = coeftest(lm.out.full.cov.incumb.comp.noec[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.comp.noec[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0 & dat.sub.comp$province2011!="Eastern Cape"),]$cat_b))
  out.mod.full.cov.opp.noec[[i]] = coeftest(lm.out.full.cov.opp.noec[[i]],vcov = vcovCluster(lm.out.full.cov.opp.noec[[i]], dat.sub[which(dat.sub$opp_control==1 & dat.sub$province2011!="Eastern Cape"),]$cat_b))
  out.mod.full.cov.opp.comp.noec[[i]] = coeftest(lm.out.full.cov.opp.comp.noec[[i]],vcov = vcovCluster(lm.out.full.cov.opp.comp.noec[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1 & dat.sub.comp$province2011!="Eastern Cape"),]$cat_b))
  
  i = i+1
}

###Generate appropriate SEs for the tables
cl.se.full.nocov = lapply(out.mod.full.nocov, `[`,,2)
cl.se.full.nocov.comp = lapply(out.mod.full.nocov.comp, `[`,,2)
cl.se.full.cov = lapply(out.mod.full.cov, `[`,,2)
cl.se.full.cov.comp = lapply(out.mod.full.cov.comp, `[`,,2)
cl.se.full.interact.cov = lapply(out.mod.full.interact.cov, `[`,,2)
cl.se.full.interact.cov.comp = lapply(out.mod.full.interact.cov.comp, `[`,,2)
cl.se.full.nocov.incumb = lapply(out.mod.full.nocov.incumb, `[`,,2)
cl.se.full.cov.incumb = lapply(out.mod.full.cov.incumb, `[`,,2)
cl.se.full.nocov.incumb.comp = lapply(out.mod.full.nocov.incumb.comp, `[`,,2)
cl.se.full.cov.incumb.comp = lapply(out.mod.full.cov.incumb.comp, `[`,,2)
cl.se.full.cov.opp = lapply(out.mod.full.cov.opp, `[`,,2)
cl.se.full.nocov.opp = lapply(out.mod.full.nocov.opp, `[`,,2)
cl.se.full.cov.opp.comp = lapply(out.mod.full.cov.opp.comp, `[`,,2)
cl.se.full.nocov.opp.comp = lapply(out.mod.full.nocov.opp.comp, `[`,,2)

cl.se.full.cov.comp.nowc = lapply(out.mod.full.cov.comp.nowc, `[`,,2)
cl.se.full.cov.comp.nokzn = lapply(out.mod.full.cov.comp.nokzn, `[`,,2)
cl.se.full.cov.comp.noec = lapply(out.mod.full.cov.comp.noec, `[`,,2)

cl.se.full.cov.nowc = lapply(out.mod.full.cov.nowc, `[`,,2)
cl.se.full.cov.nokzn = lapply(out.mod.full.cov.nokzn, `[`,,2)
cl.se.full.cov.noec = lapply(out.mod.full.cov.noec, `[`,,2)

###Table 2 in the paper: Main Effect of Change in SD on ANC Voteshare, NPE
filename = paste("npecovariatesmain.tex")
stargazer(lm.out.full.cov[[1]], lm.out.full.cov.comp[[1]],
          lm.out.full.cov[[2]], lm.out.full.cov.comp[[2]],
          lm.out.full.cov[[3]], lm.out.full.cov.comp[[3]],
          lm.out.full.cov[[4]], lm.out.full.cov.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov[[1]], cl.se.full.cov.comp[[1]],
                    cl.se.full.cov[[2]], cl.se.full.cov.comp[[2]],
                    cl.se.full.cov[[3]], cl.se.full.cov.comp[[3]],
                    cl.se.full.cov[[4]], cl.se.full.cov.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-main"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effect of change in service delivery on change on national ANC vote share",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny", 
          add.lines = list(c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),   
          out = paste(stgz,filename,sep="\\"))

###Table 3 in the paper (panel combination of the following two stargazer tables): 
#Heterogeneous Effects of Change in SD on ANC Voteshare by ANC incumbency status, NPE
filename = paste("npecovariatesANC.tex")
stargazer(lm.out.full.cov.incumb[[1]], lm.out.full.cov.incumb.comp[[1]],
          lm.out.full.cov.incumb[[2]], lm.out.full.cov.incumb.comp[[2]],
          lm.out.full.cov.incumb[[3]], lm.out.full.cov.incumb.comp[[3]],
          lm.out.full.cov.incumb[[4]], lm.out.full.cov.incumb.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.incumb[[1]], cl.se.full.cov.incumb.comp[[1]],
                    cl.se.full.cov.incumb[[2]], cl.se.full.cov.incumb.comp[[2]],
                    cl.se.full.cov.incumb[[3]], cl.se.full.cov.incumb.comp[[3]],
                    cl.se.full.cov.incumb[[4]], cl.se.full.cov.incumb.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-anconly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: ANC municipalities only",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("ANC Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

filename = paste("npecovariatesOPP.tex")
stargazer(lm.out.full.cov.opp[[1]], lm.out.full.cov.opp.comp[[1]],
          lm.out.full.cov.opp[[2]], lm.out.full.cov.opp.comp[[2]],
          lm.out.full.cov.opp[[3]], lm.out.full.cov.opp.comp[[3]],
          lm.out.full.cov.opp[[4]], lm.out.full.cov.opp.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.opp[[1]], cl.se.full.cov.opp.comp[[1]],
                    cl.se.full.cov.opp[[2]], cl.se.full.cov.opp.comp[[2]],
                    cl.se.full.cov.opp[[3]], cl.se.full.cov.opp.comp[[3]],
                    cl.se.full.cov.opp[[4]], cl.se.full.cov.opp.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-opponly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: Opposition municipalities only",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Opp Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

###APPENDIX TABLES###
###Table 15 in the Appendix:
filename = paste("npenocovariatesmain.tex")
stargazer(lm.out.full.nocov[[1]], lm.out.full.nocov.comp[[1]],
          lm.out.full.nocov[[2]], lm.out.full.nocov.comp[[2]],
          lm.out.full.nocov[[3]], lm.out.full.nocov.comp[[3]],
          lm.out.full.nocov[[4]], lm.out.full.nocov.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.nocov[[1]], cl.se.full.nocov.comp[[1]],
                    cl.se.full.nocov[[2]], cl.se.full.nocov.comp[[2]],
                    cl.se.full.nocov[[3]], cl.se.full.nocov.comp[[3]],
                    cl.se.full.nocov[[4]], cl.se.full.nocov.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-no-covariates-main"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effect of change in service delivery on change on national ANC vote share, No Covariates",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny", 
          add.lines = list(c("Covariates?","N","N","N","N","N","N","N","N"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),   
          out = paste(stgz,filename,sep="\\"))

###Table 16 in the Appendix:
filename = paste("npenocovariatesANC.tex")
stargazer(lm.out.full.nocov.incumb[[1]], lm.out.full.nocov.incumb.comp[[1]],
          lm.out.full.nocov.incumb[[2]], lm.out.full.nocov.incumb.comp[[2]],
          lm.out.full.nocov.incumb[[3]], lm.out.full.nocov.incumb.comp[[3]],
          lm.out.full.nocov.incumb[[4]], lm.out.full.nocov.incumb.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.nocov.incumb[[1]], cl.se.full.nocov.incumb.comp[[1]],
                    cl.se.full.nocov.incumb[[2]], cl.se.full.nocov.incumb.comp[[2]],
                    cl.se.full.nocov.incumb[[3]], cl.se.full.nocov.incumb.comp[[3]],
                    cl.se.full.nocov.incumb[[4]], cl.se.full.nocov.incumb.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-nocovariates-anconly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: ANC municipalities only, No Covariates",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("Covariates?","N","N","N","N","N","N","N","N"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("ANC Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

###Table 17 in the Appendix:
filename = paste("npenocovariatesOPP.tex")
stargazer(lm.out.full.nocov.opp[[1]], lm.out.full.nocov.opp.comp[[1]],
          lm.out.full.nocov.opp[[2]], lm.out.full.nocov.opp.comp[[2]],
          lm.out.full.nocov.opp[[3]], lm.out.full.nocov.opp.comp[[3]],
          lm.out.full.nocov.opp[[4]], lm.out.full.nocov.opp.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.nocov.opp[[1]], cl.se.full.nocov.opp.comp[[1]],
                    cl.se.full.nocov.opp[[2]], cl.se.full.nocov.opp.comp[[2]],
                    cl.se.full.nocov.opp[[3]], cl.se.full.nocov.opp.comp[[3]],
                    cl.se.full.nocov.opp[[4]], cl.se.full.nocov.opp.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-nocovariates-opponly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: Opposition municipalities only, No Covariates",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("Covariates?","N","N","N","N","N","N","N","N"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Opp Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

###Table 20 in the Appendix:
filename = paste("npecovariatesinteracted.tex")
stargazer(lm.out.full.interact.cov[[1]], lm.out.full.interact.cov.comp[[1]],
          lm.out.full.interact.cov[[2]], lm.out.full.interact.cov.comp[[2]],
          lm.out.full.interact.cov[[3]], lm.out.full.interact.cov.comp[[3]],
          lm.out.full.interact.cov[[4]], lm.out.full.interact.cov.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.interact.cov[[1]], cl.se.full.interact.cov.comp[[1]],
                    cl.se.full.interact.cov[[2]], cl.se.full.interact.cov.comp[[2]],
                    cl.se.full.interact.cov[[3]], cl.se.full.interact.cov.comp[[3]],
                    cl.se.full.interact.cov[[4]], cl.se.full.interact.cov.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-interacted"),
          column.separate = c(8),
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote share by local incumbent",
          #covariate.labels = covars,
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

###Table 24 in the Appendix:
filename = paste("npecovariatescompnowc.tex")
stargazer(lm.out.full.cov.comp.nowc[[1]], lm.out.full.cov.nowc[[1]],
          lm.out.full.cov.comp.nowc[[2]], lm.out.full.cov.nowc[[2]],
          lm.out.full.cov.comp.nowc[[3]], lm.out.full.cov.nowc[[3]],
          lm.out.full.cov.comp.nowc[[4]], lm.out.full.cov.nowc[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.comp.nowc[[1]], cl.se.full.cov.nowc[[1]],
                    cl.se.full.cov.comp.nowc[[2]], cl.se.full.cov.nowc[[2]],
                    cl.se.full.cov.comp.nowc[[3]], cl.se.full.cov.nowc[[3]],
                    cl.se.full.cov.comp.nowc[[4]], cl.se.full.cov.nowc[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-nowc"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effects of change in service delivery on change on national ANC vote, Excluding Western Cape",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("Covariates?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","Y", "", "Y","","Y", "", "Y", "")),
          out = paste(stgz,filename,sep="\\"))

###Table 26 in the Appendix:
filename = paste("npecovariatescompnokzn.tex")
stargazer(lm.out.full.cov.comp.nokzn[[1]], lm.out.full.cov.nokzn[[1]],
          lm.out.full.cov.comp.nokzn[[2]], lm.out.full.cov.nokzn[[2]],
          lm.out.full.cov.comp.nokzn[[3]], lm.out.full.cov.nokzn[[3]],
          lm.out.full.cov.comp.nokzn[[4]], lm.out.full.cov.nokzn[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.comp.nokzn[[1]], cl.se.full.cov.nokzn[[1]],
                    cl.se.full.cov.comp.nokzn[[2]], cl.se.full.cov.nokzn[[2]],
                    cl.se.full.cov.comp.nokzn[[3]], cl.se.full.cov.nokzn[[3]],
                    cl.se.full.cov.comp.nokzn[[4]], cl.se.full.cov.nokzn[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-nokzn"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effects of change in service delivery on change on national ANC vote, Excluding KwaZulu-Natal",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(c("Covariates?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","Y", "", "Y","","Y", "", "Y", "")),
          out = paste(stgz,filename,sep="\\"))

###Table 28 in the Appendix:
filename = paste("npecovariatescompnoec.tex")
stargazer(lm.out.full.cov.comp.noec[[1]], lm.out.full.cov.noec[[1]],
          lm.out.full.cov.comp.noec[[2]], lm.out.full.cov.noec[[2]],
          lm.out.full.cov.comp.noec[[3]], lm.out.full.cov.noec[[3]],
          lm.out.full.cov.comp.noec[[4]], lm.out.full.cov.noec[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.comp.noec[[1]], cl.se.full.cov.noec[[1]],
                    cl.se.full.cov.comp.noec[[2]], cl.se.full.cov.noec[[2]],
                    cl.se.full.cov.comp.noec[[3]], cl.se.full.cov.noec[[3]],
                    cl.se.full.cov.comp.noec[[4]], cl.se.full.cov.noec[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-noec"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effects of change in service delivery on change on national ANC vote, Excluding Easternn Cape",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny", 
          add.lines = list(c("Covariates?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("$<$80th Percentile SD?","Y","Y","Y","Y","Y","Y","Y","Y"),
                           c("Competitive only?","Y", "", "Y","","Y", "", "Y", "")),
          out = paste(stgz,filename,sep="\\"))

###################################
######INCLUDING TOP QUINTILE#######
###################################
i = 1
out.mod.full.cov = list()
out.mod.full.cov.comp = list()
out.mod.full.cov.incumb = list()
out.mod.full.cov.incumb.comp = list()
out.mod.full.cov.opp = list()
out.mod.full.cov.opp.comp = list()

lm.out.full.cov = list()
lm.out.full.cov.comp = list()
lm.out.full.cov.incumb = list()
lm.out.full.cov.incumb.comp = list()
lm.out.full.cov.opp = list()
lm.out.full.cov.opp.comp = list()

for(iv in c("water_ch10", "flush_ch10", "refuse_ch10", "sd_means_ch10")){
  dat.sub = dat[,c(iv,
                   "npeanc_vs_ch10", "unempfrac_ch10", "log_income_ch10", "log_pop_ch10", "femalefrac_ch10", "blackfrac_ch10", 
                   "whitefrac_ch10", "indianfrac_ch10", "colouredfrac_ch10", "traddwell_ch10", "cat_b", "opp_control", "province2011",
                   "employment_broad2001","log_income2001","log_pop2001","sex2001","indian_frac2001","colored_frac2001","black_frac2001","white_frac2001",
                   "share_area_all_tbvc_ta","traddwell_proportion2001","refuse_perc2001","all_flush_perc2001","home_water_perc2001","sd_means2001",
                   "competitive")]
  dat.sub = na.omit(dat.sub)
  
  dat.sub.comp = dat.sub[which(dat.sub$competitive==1),]
  
  dat.sub.comp = na.omit(dat.sub.comp)
  
  iv.int = paste(iv,"*opp_control", sep="")
  
  if(iv=="water_ch10"){
    baseline = " + home_water_perc2001"
  } else if(iv=="elec_ch10"){
    baseline = " + elec_perc2001"
  } else if(iv=="flush_ch10"){
    baseline = " + all_flush_perc2001"
  } else if(iv=="refuse_ch10"){
    baseline = " + refuse_perc2001"
  } else {
    baseline = " + sd_means2001"
  }
  
  #######REGRESSIONS:
  ##Main specifications
  lm.out.full.cov[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, cov, baseline, sep=""), data=dat.sub)
  lm.out.full.cov.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv, cov, baseline, sep=""), data=dat.sub.comp)
  
  ##By incumbency status
  lm.out.full.cov.incumb[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==0),])
  lm.out.full.cov.incumb.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==0),])
  
  lm.out.full.cov.opp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub[which(dat.sub$opp_control==1),])
  lm.out.full.cov.opp.comp[[i]] = lm(paste("npeanc_vs_ch10", " ~ ", iv.int, cov, baseline, sep=""), data=dat.sub.comp[which(dat.sub.comp$opp_control==1),])
    
  #######CLUSTER STANDARD ERRORS:
  out.mod.full.cov[[i]] = coeftest(lm.out.full.cov[[i]],vcov = vcovCluster(lm.out.full.cov[[i]], dat.sub$cat_b))
  out.mod.full.cov.comp[[i]] = coeftest(lm.out.full.cov.comp[[i]],vcov = vcovCluster(lm.out.full.cov.comp[[i]], dat.sub.comp$cat_b))

  out.mod.full.cov.incumb[[i]] = coeftest(lm.out.full.cov.incumb[[i]],vcov = vcovCluster(lm.out.full.cov.incumb[[i]], dat.sub[which(dat.sub$opp_control==0),]$cat_b))
  out.mod.full.cov.incumb.comp[[i]] = coeftest(lm.out.full.cov.incumb.comp[[i]],vcov = vcovCluster(lm.out.full.cov.incumb.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==0),]$cat_b))
  
  out.mod.full.cov.opp[[i]] = coeftest(lm.out.full.cov.opp[[i]],vcov = vcovCluster(lm.out.full.cov.opp[[i]], dat.sub[which(dat.sub$opp_control==1),]$cat_b))
  out.mod.full.cov.opp.comp[[i]] = coeftest(lm.out.full.cov.opp.comp[[i]],vcov = vcovCluster(lm.out.full.cov.opp.comp[[i]], dat.sub.comp[which(dat.sub.comp$opp_control==1),]$cat_b))
  
  i = i+1
}

###Generate appropriate SEs for the tables
cl.se.full.cov = lapply(out.mod.full.cov, `[`,,2)
cl.se.full.cov.comp = lapply(out.mod.full.cov.comp, `[`,,2)
cl.se.full.cov.incumb = lapply(out.mod.full.cov.incumb, `[`,,2)
cl.se.full.cov.incumb.comp = lapply(out.mod.full.cov.incumb.comp, `[`,,2)
cl.se.full.cov.opp = lapply(out.mod.full.cov.opp, `[`,,2)
cl.se.full.cov.opp.comp = lapply(out.mod.full.cov.opp.comp, `[`,,2)

###Appendix table 10: Main results including top quintile, NPE
filename = paste("topqincludednpecovariatesmain.tex")
stargazer(lm.out.full.cov[[1]], lm.out.full.cov.comp[[1]],
          lm.out.full.cov[[2]], lm.out.full.cov.comp[[2]],
          lm.out.full.cov[[3]], lm.out.full.cov.comp[[3]],
          lm.out.full.cov[[4]], lm.out.full.cov.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov[[1]], cl.se.full.cov.comp[[1]],
                    cl.se.full.cov[[2]], cl.se.full.cov.comp[[2]],
                    cl.se.full.cov[[3]], cl.se.full.cov.comp[[3]],
                    cl.se.full.cov[[4]], cl.se.full.cov.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-main"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Effect of change in service delivery on change on national ANC vote share",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny", 
          add.lines = list(
            c("Competitive only?","","Y","","Y","","Y","","Y")),   
          out = paste(stgz,filename,sep="\\"))

###Appendix table 11: Main results including top quintile, ANC only, NPE
filename = paste("topqincludednpecovariatesANC.tex")
stargazer(lm.out.full.cov.incumb[[1]], lm.out.full.cov.incumb.comp[[1]],
          lm.out.full.cov.incumb[[2]], lm.out.full.cov.incumb.comp[[2]],
          lm.out.full.cov.incumb[[3]], lm.out.full.cov.incumb.comp[[3]],
          lm.out.full.cov.incumb[[4]], lm.out.full.cov.incumb.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.incumb[[1]], cl.se.full.cov.incumb.comp[[1]],
                    cl.se.full.cov.incumb[[2]], cl.se.full.cov.incumb.comp[[2]],
                    cl.se.full.cov.incumb[[3]], cl.se.full.cov.incumb.comp[[3]],
                    cl.se.full.cov.incumb[[4]], cl.se.full.cov.incumb.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-anconly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: ANC municipalities only",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(
            c("ANC Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
            c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))

###Appendix table 12: Main results including top quintile, Opp only, NPE
filename = paste("topqincludednpecovariatesOPP.tex")
stargazer(lm.out.full.cov.opp[[1]], lm.out.full.cov.opp.comp[[1]],
          lm.out.full.cov.opp[[2]], lm.out.full.cov.opp.comp[[2]],
          lm.out.full.cov.opp[[3]], lm.out.full.cov.opp.comp[[3]],
          lm.out.full.cov.opp[[4]], lm.out.full.cov.opp.comp[[4]],
          #This next row includes the correct SEs:
          se = list(cl.se.full.cov.opp[[1]], cl.se.full.cov.opp.comp[[1]],
                    cl.se.full.cov.opp[[2]], cl.se.full.cov.opp.comp[[2]],
                    cl.se.full.cov.opp[[3]], cl.se.full.cov.opp.comp[[3]],
                    cl.se.full.cov.opp[[4]], cl.se.full.cov.opp.comp[[4]]),
          style = "ajps",
          label = paste("tab:npe-covariates-opponly"),
          column.separate = c(8),
          #covariate.labels = covars,
          title = "Heterogeneous effects of change in service delivery on change on national ANC vote: Opposition municipalities only",
          model.numbers=F,
          dep.var.labels.include = FALSE,
          column.labels = c("Change in ANC Vote Share (National Elections)"),
          #column.separate = c(1,1,1,1,1),
          omit.stat = c("f", "adj.rsq", "ser"),
          notes = "Standard errors clustered by municipality",
          float.env = "sidewaystable",
          font.size = "tiny",
          add.lines = list(
            c("Opp Control Only?","Y","Y","Y","Y","Y","Y","Y","Y"),
            c("Competitive only?","","Y","","Y","","Y","","Y")),
          out = paste(stgz,filename,sep="\\"))
